For the sensitivity testing we estimated the size spectrum slope and intercept for every year and within each survey area an additional time, at each pass one of the species was omitted from the analysis. This was done to capture the relative importance each species had on the estimation of the biomass size spectra.
The biological data used as a starting point is the same “node” used for the size spectrum estimations using the full community. This is loaded so that we can interrogate the overall biomass patterns that correspond to the results we see from the sensitivity testing.
# Access Biological Data with {gmRi}
# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")
withr::with_dir(rprojroot::find_root('_targets.R'),
tar_load(nefsc_1g_binned))
# Get Yearly and Regionally averaged Biomasses
yr_szn_overall <- nefsc_1g_binned %>%
group_by(Year, survey_area) %>%
summarise(overall_bio = sum(strat_total_lwbio_s, na.rm = T),
overall_abund = sum(strat_total_abund_s, na.rm = T),
.groups = "drop")
# Get what it is for each species
sp_yr_szn <- nefsc_1g_binned %>%
group_by(comname, spec_class, fishery, hare_group, Year, survey_area) %>%
summarise(sp_bio = sum(strat_total_lwbio_s, na.rm = T),
sp_abund = sum(strat_total_abund_s, na.rm = T),
.groups = "drop")
# Join in the overall, get the fraction and what the percent would be
perc_bio <- left_join(sp_yr_szn, yr_szn_overall, by = c("Year", "survey_area")) %>%
mutate(
bio_frac = sp_bio / overall_bio,
abund_frac = sp_abund / overall_abund,
bio_perc = bio_frac * 100,
abund_perc = abund_frac * 100
)
# perc_bio overall rankings
region_overalls <- perc_bio %>%
group_by(comname, spec_class, fishery, hare_group, survey_area) %>%
summarise(avg_bio_perc = mean(bio_perc, na.rm = T),
avg_abund_perc = mean(abund_perc, na.rm = T),
.groups = "drop")
# Treemaps
ggplot(region_overalls) +
geom_treemap(aes(area = avg_bio_perc, fill = comname)) +
geom_treemap_text(aes(area = avg_bio_perc, fill = comname, label = comname)) +
facet_wrap(~survey_area, ncol = 2) +
scale_fill_gmri() +
theme(legend.position = "none") +
labs(title = "Fraction of Stratified Biomass: 1970-2019")
# Treemaps
ggplot(region_overalls) +
geom_treemap(aes(area = avg_abund_perc, fill = comname)) +
geom_treemap_text(aes(area = avg_abund_perc, fill = comname, label = comname)) +
facet_wrap(~survey_area, ncol = 2) +
scale_fill_gmri() +
theme(legend.position = "none") +
labs(title = "Fraction of Stratified Abundance: 1970-2019")
# Play around with plotting yearly influences
ggplot(perc_bio, aes(Year, bio_perc)) +
# geom_point(aes(color = comname), show.legend = F) +
# geom_line(aes(color = comname, group = comname), show.legend = F) +
geom_col(aes(fill = spec_class)) +
facet_grid(survey_area~.) +
scale_color_gmri() +
scale_fill_gmri() +
theme(legend.position = "bottom") +
labs(y = "Percent of Yearly Biomass",
x = NULL,
fill = "")
#
# Play around with plotting yearly influences
ggplot(perc_bio, aes(Year, abund_perc)) +
# geom_point(aes(color = comname), show.legend = F) +
# geom_line(aes(color = comname, group = comname), show.legend = F) +
geom_col(aes(fill = spec_class)) +
facet_grid(survey_area~.) +
#facet_grid(spec_class~survey_area) +
scale_color_gmri() +
scale_fill_gmri() +
theme(legend.position = "bottom") +
labs(y = "Percent of Yearly Biomass",
x = NULL,
fill = "")
#
# Access Biological Data with {gmRi}
# nmfs_bio <- gmRi::gmri_survdat_prep(survdat_source = "bio")
withr::with_dir(rprojroot::find_root('_targets.R'),
tar_load(species_sensitivity_shifts))
# Make year a numeric value
sens_results <- species_sensitivity_shifts %>%
mutate(Year = as.numeric(Year))
# Join in the biomass info
sens_results <- sens_results %>%
select(-spec_class) %>%
left_join(perc_bio,
by = c("Year" = "Year",
"survey_area" = "survey_area",
"spec_omit" = "comname"))
sens_ranks <- sens_results %>%
group_by(
`Species Removed` = spec_omit,
Area = survey_area,
Guild = spec_class,
Fishery = fishery) %>%
summarise(`Avg. Slope Shift` = mean(slope_shift, na.rm = T),
`Avg. Slope Shift (z)` = mean(slope_shift_z, na.rm = T),
.groups = "drop") %>%
mutate(across(where(is.numeric), round, 2))
# Slope Steepeners
sens_ranks %>%
arrange(`Avg. Slope Shift (z)`) %>%
head(8) %>%
gt::gt() %>%
gt::tab_header(
title = "Top 8 Steepeners",
subtitle = "Removing these species causes slope to steepen"
)
| Top 8 Steepeners | |||||
|---|---|---|---|---|---|
| Removing these species causes slope to steepen | |||||
| Species Removed | Area | Guild | Fishery | Avg. Slope Shift | Avg. Slope Shift (z) |
| bullnose ray | SNE | NA | Not Labelled | -0.41 | -2.16 |
| silver hake | GoM | Groundfish | Commercially Targeted | -0.20 | -1.08 |
| silver hake | SNE | Groundfish | Commercially Targeted | -0.14 | -0.75 |
| barndoor skate | GB | Elasmobranch | Commercially Targeted | -0.11 | -0.59 |
| butterfish | MAB | Pelagic | Not Commercially Targeted | -0.08 | -0.58 |
| silver hake | GB | Groundfish | Commercially Targeted | -0.11 | -0.56 |
| roughtail stingray | SNE | NA | Not Labelled | -0.09 | -0.49 |
| spotted hake | MAB | Groundfish | Not Commercially Targeted | -0.06 | -0.44 |
# Slope Flatteners
sens_ranks %>%
arrange(desc(`Avg. Slope Shift (z)`)) %>%
head(8) %>%
gt::gt() %>%
gt::tab_header(
title = "Top 8 Flatteners",
subtitle = "Removing these species causes slope to flatten"
)
| Top 8 Flatteners | |||||
|---|---|---|---|---|---|
| Removing these species causes slope to flatten | |||||
| Species Removed | Area | Guild | Fishery | Avg. Slope Shift | Avg. Slope Shift (z) |
| spiny dogfish | MAB | Elasmobranch | Commercially Targeted | 0.13 | 0.92 |
| spiny dogfish | SNE | Elasmobranch | Commercially Targeted | 0.15 | 0.80 |
| winter skate | GB | Elasmobranch | Commercially Targeted | 0.04 | 0.22 |
| spiny dogfish | GoM | Elasmobranch | Commercially Targeted | 0.04 | 0.21 |
| atlantic cod | GB | Groundfish | Commercially Targeted | 0.04 | 0.19 |
| spiny dogfish | GB | Elasmobranch | Commercially Targeted | 0.03 | 0.17 |
| atlantic cod | GoM | Groundfish | Commercially Targeted | 0.02 | 0.12 |
| sand tiger | MAB | NA | Not Labelled | 0.02 | 0.12 |
# Do we need to show the relative influence?
# Things seem to wash out with the raw z-scores
# Maybe they should be adjusted/weighted for how much any single omission
# pulls a given year compared to either the
# average z-shift for that year+region
# or by the max that it shifts....
# plot heatmap
influence_heatmap <- function(sens_results,
filt_area,
lims = c(-1, 1),
rescale_z = F){
# Make breaks to indicate limits are truncating
blabs <- seq(from = lims[1], to = lims[2], length.out = 5)
blabs[c(1,5)] <- c(str_c(">", blabs[1]),
str_c("<", blabs[5]))
# Filter data
filt_dat <- sens_results %>%
mutate(Year = as.numeric(Year)) %>%
filter(survey_area == filt_area)
# Rescale?
if(rescale_z == T){
message("Rescaling not finished")
}
# Make heatmap
filt_dat %>%
ggplot(aes(Year, fct_rev(spec_omit), fill = slope_shift_z)) +
geom_tile() +
scale_fill_distiller(palette = "RdBu",
limits = lims,
labels = blabs,
oob = scales::oob_squish,
na.value = "gray20") +
scale_x_continuous(expand = expansion(add = c(0,0))) +
facet_wrap(~survey_area, ncol = 2, scales = "free") +
#facet_grid(spec_class~survey_area, scales = "free") +
#facet_grid(fishery~survey_area, scales = "free") +
theme(legend.position = "bottom",
legend.title = element_text(angle = 0),
axis.text.y = element_text(size = 7),
panel.background = element_rect(fill = "gray40"),
panel.grid = element_line(color = "transparent")) +
guides(fill = guide_colorbar(title.position = "top",
title.hjust = 0.5,
barwidth = unit(3.25, "in"))) +
labs(y = "Species omitted",
fill = "Shift in Slope (z-score)")
}
# GOM
sens_results %>%
#filter(spec_class == "Groundfish") %>%
influence_heatmap(filt_area = "GoM", lims = c(-2,2))
# GOM
influence_heatmap(sens_results, "GB", lims = c(-2,2))
# GOM
influence_heatmap(sens_results, "SNE", lims = c(-4,4))
# GOM
influence_heatmap(sens_results, "MAB", lims = c(-4,4))
sens_results %>%
ggplot(aes(bio_frac, slope_shift_z)) +
geom_point(alpha = 0.3) +
geom_smooth(formula = "y ~ x", method = "lm") +
labs(x = "Fration of Biomass", y = "Shift in Spectra Slope (z)")
I think this might help get around the “wash out” of the heatmap. The color scale should be longer, or potentially just not center on white…
# plot a species:
# multiple panel plot function
# Function for plotting a species that stands out
plot_spectra_influence <- function(sens_results, filter_spec, filter_area){
# Grab the relevant data
filt_d <- sens_results %>%
filter(spec_omit == filter_spec,
survey_area == filter_area)
# Make Plot Label
p_label <- str_c(str_to_title(filter_spec), " Omission")
# 1. Actual Value Shift
# 1a. Slope
splot <- filt_d %>%
ggplot(aes(Year, l10_slope_strat)) +
geom_line(aes(color = "All Species"), group = 1) +
geom_point(aes(color = "All Species")) +
geom_line(aes(y = omit_slope_strat, color = p_label), group = 1) +
geom_point(aes(y = omit_slope_strat, color = p_label)) +
facet_wrap(~survey_area, ncol = 1) +
scale_color_gmri() +
theme(legend.position = "bottom") +
labs(color = NULL,
y = "Biomass Spectrum Slope",
x = NULL)
# 2. What the z-score shift is
z_shift <- filt_d %>%
ggplot(aes(Year, slope_shift_z)) +
geom_line() +
geom_point() +
labs(y = "Slope Shift (z)",
x = NULL)
# Biomass Fraction Each Year:
bfrac <- filt_d %>%
ggplot(aes(Year, bio_frac)) +
geom_line() +
geom_point() +
scale_y_continuous(labels = scales::percent) +
labs(y = "Biomass Fraction")
# Do they fit
plotsemble <- splot / z_shift / bfrac
return(plotsemble)
}
plot_spectra_influence(sens_results = sens_results,
filter_spec = "haddock",
filter_area = "GoM")
plot_spectra_influence(sens_results = sens_results,
filter_spec = "atlantic herring",
filter_area = "SNE")
plot_spectra_influence(sens_results = sens_results,
filter_spec = "silver hake",
filter_area = "GoM")
plot_spectra_influence(sens_results = sens_results,
filter_spec = "silver hake",
filter_area = "GB")
plot_spectra_influence(sens_results = sens_results,
filter_spec = "butterfish",
filter_area = "GoM")
plot_spectra_influence(sens_results = sens_results,
filter_spec = "atlantic cod",
filter_area = "GoM")